// Simulates the effect of implementing methane pricing
// Author Levi Marks levi.marks.1>at>gmail.com

// 0 // Preliminaries
version 13 // Version control for random number generation
global DATA "[Insert Path Here]"
global iterations 1000
set more off
graph drop _all

// 1 // Setup, trimming dataset

// 1.1 // Create compacted dataset to speed simulation

use "$DATA/ghgrp_di_intermediate", clear

order rate p_spot
sort p_spot facility_id year rate 

keep if trim_5 != 1

gen p_fp1 = p_spot^(-2)
gen p_fp2 = p_spot^(-2)*log(p_spot)

keep rate p_spot p_fp1 p_fp2 facility_id year region_year wells oil gas completions parent 
sort facility_id year
bysort facility_id: gen group_id = _n == 1
replace group_id = sum(group_id)
	
compress*
save "$DATA/Temp/full_compressed.dta", replace
*/

// 1.2 // Initial regression to recover parameters (needed for most sections other than combining results)

use "$DATA/ghgrp_di_intermediate", clear

// Preliminary regression macros
local dependent			rate
local independent  		p_spot 
local controls		 	wells oil completions
local FE 				i.region_year i.facility_id 
local conditions		if trim_5 != 1

// Second-order FP Regression
fp <`independent'>, dimension(2) powers(-2 -2) replace: reg `dependent' <`independent'> `controls' `FE' `weights' `conditions', vce(cluster parent)

// Recovering parameters
mat coef = e(b)
mat fp_pwr =  e(fp_fp)
scalar beta_1 = coef[1,1]
scalar beta_2 = coef[1,2]
scalar constant = coef[1,`=colsof(coef)']
scalar fp_pwr_1 = fp_pwr[1,1]
scalar fp_pwr_2 = fp_pwr[1,2]
*/



// 2 // Graphs

// 2.1 // Setup for graphs

// More parameters to adjust graph
qui sum rate `conditions'
local avg_rate `=r(mean)'
qui sum p_spot `conditions'
local avg_price `=r(mean)'
local rate_midpoint = coef[1,`=colsof(coef)'] + coef[1,1]*(`avg_price'^`=fp_pwr[1,1]') + coef[1,2]*log(`avg_price')*(`avg_price'^`=fp_pwr[1,1]')
local adjust `=`rate_midpoint' - `avg_rate'' // This is used to adjust the graph y-axis to center the regression line at the average emission rate

keep `conditions'
sum p_spot `conditions'
local mean_p = r(mean)
local min_p = r(min)
local max_p = r(max)

collapse (mean) rate gas p_spot, by(facility_id)

egen total_gas = sum(gas) // Total gas in sample
replace gas = gas*(32592000000/total_gas) // Scaling up gas production based on EIA 2016 estimate (from EIA March 2020 Monthly Energy Report)
gen emissions = rate*gas
sum emissions
scalar total_emissions = r(sum)

egen lower_bound = min(rate)
order lower_bound, after(facility_id)
rename rate rate_0
rename gas gas_0
rename p_spot price_0
gen emissions_0 = rate_0*gas_0
gen ag_cost_0 = 0
gen ag_reduction_0 = 0
gen ag_value_0 = 0
gen rate_00 = rate_0 // Reference so others can be deleted
gen price_00 = price_0 // Reference so others can be deleted
*/

// Note: Do not comment out the code above this line

// 2.2 // Graph to illustrate the simulation

// Simulation emission reductions
local iterations 224
forval i = 1/`iterations' {
	di `i'
	gen price_`i' = price_0 + .79*.05*`i'
	gen rate_`i'= rate_`=(`i'-1)' + .79*.05*(`=fp_pwr_1'*beta_1*price_`i'^`=fp_pwr_1-1' + beta_2*( `=fp_pwr_2'*(price_`i'^`=fp_pwr_2-1')*log(price_`i') + (price_`i'^`=fp_pwr_2')*(1/price_`i')  ) ) // First derivative
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_0 if rate_0 < lower_bound // Reference start point
}

// Constructing graph
keep facility_id rate_* price_*
reshape long rate_ price_, i(facility_id) j(step)
rename rate_ rate
rename price_ price

// For arrows
bysort facility_id: gen rate_y1 = rate[`iterations'] // Adding 1 because there is a zero
bysort facility_id: gen price_x1 = price[`iterations']
bysort facility_id: gen rate_y2 = rate[`=`iterations'+1']
bysort facility_id: gen price_x2 = price[`=`iterations'+1']

levelsof facility_id, local(facilities)
foreach facility in `facilities' {
	local lines `lines' line rate price if facility_id == `facility' & rate < .02, lcolor(gs10) lwidth(vthin) || 
}

twoway line rate price if facility_id == 0 || ///
	`lines' ///
	function y = constant - `adjust' + beta_1*x^(fp_pwr_1) + beta_2*x^(fp_pwr_2)*log(x), range(1.25 `=3.87+.79*11.2') lwidth(medthick) lcolor(gs3) lpattern(dash) || ///
	function y = constant - `adjust' + beta_1*x^(fp_pwr_1) + beta_2*x^(fp_pwr_2)*log(x), range(1.25 5.79) lwidth(medthick) lcolor(gs3) || ///
	scatter rate price if rate < .02 & step == 0, msymbol(plus) msize(medium) mcolor(gs8) || ///
	pcarrow rate_y1 price_x1 rate_y2 price_x2 if step == 0 & rate_y2 < .02, mcolor(gs10) lcolor(gs10) mlwidth(vthin) msize(1.5) barbsize(.25) lwidth(vthin) ///
	graphregion(color(white)) xsize(12) ysize(4) ///
	xlabel(0(2)12, labsize(medlarge) ) xtitle("Price ($/Mcf)", size(large) margin(0 0 -2 2)) ///
	ytitle("Emission Rate", size(large) margin(0 2 0 0)) ///
	yscale(range(-.0001 .0201)) ylabel(0(.005).02, labsize(medlarge) glcolor(gs14) glwidth(vvthin) ) ///
	legend(order(216 "Fitted" "Relationship" 215 "Continuation" "past Support" "for Prices" 217 "Facility" "Start Point" 218 "Predicted" "Rate for $20" "CO2 Price" ) ///
	size(medlarge) region(lwidth(thin) lcolor(gs14) margin(1 1 1 1)) symxsize(5) margin(1 1 1 1) bmargin(5 1 1 1) pos(2) col(1))
graph export "$DATA/Results/p8_sim.pdf", replace
a
*/

// 2.3 // MACC (comment out 2.2 before running this)
/*
local iterations 670
local step_size .05

forval i = 1/`iterations' {
	di `i'
	gen price_`i' = price_00 + .79*`step_size'*`i'
	gen rate_`i'= rate_`=(`i'-1)' + .79*`step_size'*(`=fp_pwr_1'*beta_1*price_`i'^`=fp_pwr_1-1' + beta_2*( `=fp_pwr_2'*(price_`i'^`=fp_pwr_2-1')*log(price_`i') + (price_`i'^`=fp_pwr_2')*(1/price_`i')  ) ) // First derivative
	qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
	qui replace rate_`i' = rate_00 if rate_00 < lower_bound // Reference start point
	
	gen emissions_`i' = rate_`i'*gas_0
	gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
	gen value_`i' = reduction_`i'*price_00 // Value of conserved gas
	gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
	
	egen sum_cost_`i' = sum(cost_`i')
	egen sum_reduction_`i' = sum(reduction_`i')
	egen sum_value_`i' = sum(value_`i')
	gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
	gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
	gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
	gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
	
	drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
}

keep ag_* mc_*
keep if _n == 1

gen dummy = 1
reshape long ag_cost_ ag_reduction_ ag_value_ mc_, i(dummy) j(tax)
drop dummy

rename tax tax_mcf
rename mc_ mc_mcf
rename ag_reduction_ ag_reduction_mcf
rename ag_cost_ ag_cost
rename ag_value_ ag_value
gen marginal_value = 3.123 // Calculated this by setting the step size to .001. Using this for the graph to create a straight line for readability; however, point estimates in the paper are based on the "ag_value" variable, which will vary slightly from the average depending on which facilities reach the lower bound first.

replace mc_mcf = marginal_value if tax_mcf == 0

replace tax_mcf = tax_mcf*`step_size'
recast double tax_mcf
replace tax_mcf = round(tax_mcf,.01)
gen tax_tco2 = (tax*53.68)/34
gen ag_reduction_tco2e = ag_reduction_mcf*34/53.68
gen mc_tco2e = mc_mcf*53.68/34
gen reduction_percent = (ag_reduction_mcf/total_emissions)*100
gen net_cost = (ag_cost - ag_value)/28400000000

replace ag_reduction_mcf = 0 if missing(ag_reduction_mcf)
replace ag_reduction_tco2e = 0 if missing(ag_reduction_tco2e)

format mc* tax* reduction %6.2fc
format ag* %12.0fc

order tax_mcf tax_tco2 mc_mcf mc_tco2e ag_cost ag_reduction_mcf ag_reduction_tco2e reduction_percent

// Interpolating $5, $20, and $44.4 carbon prices
set obs `=_N+1'
replace tax_tco2 = 5 if missing(tax_tco2)
set obs `=_N+1'
replace tax_tco2 = 20 if missing(tax_tco2)
set obs `=_N+1'
replace tax_tco2 = 44.4 if missing(tax_tco2)
sort tax_tco2
gen percent_next = (tax_tco2 - tax_tco2[_n-1])/(tax_tco2[_n+1] - tax_tco2[_n-1]) if missing(mc_mcf)
gen percent_prev = (tax_tco2[_n+1] - tax_tco2)/(tax_tco2[_n+1] - tax_tco2[_n-1]) if missing(mc_mcf)
replace tax_mcf = tax_mcf[_n-1]*percent_prev + tax_mcf[_n+1]*percent_next if missing(tax_mcf)
replace tax_tco2 = tax_tco2[_n-1]*percent_prev + tax_tco2[_n+1]*percent_next if missing(tax_tco2)
foreach var in mc_mcf mc_tco2e ag_cost ag_reduction_mcf ag_reduction_tco2e reduction_percent ag_value net_cost {
	replace `var' = `var'[_n-1]*percent_prev + `var'[_n+1]*percent_next if missing(`var')
}

// Value of conserved gas for area plot
set obs `=_N+1'
replace marginal_value = marginal_value[_n-1] if missing(marginal_value)
replace reduction_percent = 100 if missing(reduction_percent)

// Indicating out-of-sample predictions: Only 79% that is methane is taxed, so the highest observe price (5.75) is equivalent to a methane tax of 5.75/.79 = 7.278
set obs `=_N+1'
replace reduction_percent = 97.75 if _n == _N
gen reduction_percent_end = 97.75 if _n == _N
replace mc_mcf = 7.278 if _n == _N
gen mc_mcf_end = 29.75 if _n == _N

set obs `=_N+2' // This makes a little dash
replace mc_mcf = 7.278 if _n >= _N - 1
replace reduction_percent = 97 if _n == _N - 1
replace reduction_percent = 98.5 if _n == _N

twoway area marginal_value reduction_percent, color(gs14) lwidth(vvthin) lcolor(gs14) || ///
	line mc_mcf reduction_percent if _n < _N - 3, lcolor(black) lwidth(vthin) || /// Thin line then thick line to double check axes are scaled correctly
	line mc_tco2e ag_reduction_tco2e, lcolor(gs3) lwidth(medthick) yaxis(2) xaxis(2) || ///
	line mc_mcf reduction_percent if _n >= _N - 1, lcolor(gs12) lwidth(thin) || ///
	pcarrow mc_mcf reduction_percent mc_mcf_end reduction_percent_end if _n == _N - 2, lcolor(gs12) mcolor(gs12) mlwidth(thin) msize(1.75) barbsize(.35) lwidth(thin) ///
	yline(5.68, lcolor(gs8) lpattern(dash)) /// Obtain these by examining mc_mcf in the dataset at i.e. tax_tco2 = 5
	yline(12.67, lcolor(gs8) lpattern(dash)) ///
	yline(25.41, lcolor(gs8) lpattern(dash)) ///
	text(5.68 2 "$5 Carbon Price" , yaxis(1) xaxis(1) placement(ne) orient(hor) margin(.5 .5 .5 .5) color(gs8) size(medsmall) ) ///
	text(12.67 2 "$20 Carbon Price", yaxis(1) xaxis(1) placement(ne) orient(hor) margin(.5 .5 .5 .5) color(gs8) size(medsmall) ) ///
	text(25.41 2 "Social Cost of Methane", yaxis(1) xaxis(1) placement(ne) orient(hor) margin(.5 .5 .5 .5) color(gs8) size(medsmall) ) ///
	text(.5 50 "Value of Captured Gas", yaxis(1) xaxis(1) placement(n) orient(hor) margin(.5 .5 .5 .5) color(gs8) size(medsmall) ) ///
	text(`=(29.75+7.278)/2' 97.5 "Out-of-Sample" "Predictions" , yaxis(1) xaxis(1) placement(w) orient(vert) margin(.5 .5 .5 .5) color(gs12) size(*.8)) ///
	graphregion(color(white)) legend(off) xsize(6) ysize(4) ///
	plotregion(margin(0 0 0 0)) ///
	xlabel(0 "0%" 20 "20%" 40 "40%" 60 "60%" 80 "80%" 100 "100%", axis(1) labsize(medsmall) format(%12.0fc) nogrid) ///
	xlabel(15000000(15000000)60000000 , axis(2) labsize(medsmall) format(%12.0fc) nogrid labcolor(gs8) tlcolor(gs8) tlwidth(medthin)) ///
	xscale(range(0 100)) ///
	xscale(range(0 67125591) axis(2) lcolor(gs8) lwidth(medthin)) /// Calculate total CO2e emissions for this one
	xtitle("Total Abatement"  , axis(1) size(medlarge) margin(0 0 -2 2)) ///
	xtitle("Total Abatement (tCO2e)", axis(2) size(medlarge) margin(0 0 2 -2) color(gs8)) ///
	ylabel(0(5)31, axis(1) labsize(medsmall)format(%6.0fc) nogrid) ///
	yscale(range(0 31)) ///
	ylabel(10(10)40, axis(2) labsize(medsmall) format(%6.0fc) nogrid labcolor(gs8) tlcolor(gs8) tlwidth(medthin)) ///
	yscale(range(0 `=31*53.68/34') axis(2) lcolor(gs8) lwidth(medthin)) ///
	ytitle("Marginal Cost ($/Mcf)"  , axis(1) size(medlarge) margin(0 2 0 0)) ///
	ytitle("Marginal Cost ($/tCO2e)", axis(2) size(medlarge) margin(2 0 0 0) color(gs8)) 
graph export "$DATA/Results/p8_macc.pdf", replace
a
*/



// 3 // Bootstrap for coefficient estimates with SEs
/*
set seed 79
local iterations $iterations

// 3.1 // Construct matrix of bootstrapped coefficient estimates

use "$DATA/Temp/full_compressed.dta", clear

// Preliminary regression macros
local dependent			rate
local independent  		p_spot 
local controls		 	wells oil completions
local FE 				i.region_year i.facility_id 
local conditions		if trim_5 != 1

// Initial regression

areg `dependent' p_fp1 p_fp2 `controls' i.region_year, absorb(facility_id) vce(cluster parent)
mat Coef = e(b)
mat Results = Coef[1,1],Coef[1,2],Coef[1,`=colsof(Coef)']
mat colnames Results = b_1 b_2 constant

// Bootstrap regressions
forval k = 1/`iterations' {
	di `k'
	// Creating bootstrap sample
	use "$DATA/Temp/full_compressed.dta", clear
	scalar select = ceil(runiform()*213)
	
	qui mkmat * if group_id == select, matrix(Boot_Sample)
	
	forval i = 1/212 {
		scalar select = ceil(runiform()*213)
		mkmat * if group_id == select, matrix(New)
		mat Boot_Sample = Boot_Sample\New
	}
	
	drop *
	qui svmat Boot_Sample, names(col)
	
	// Regression and saving results
	qui areg `dependent' p_fp1 p_fp2 `controls' i.region_year, absorb(facility_id) vce(cluster parent)
	mat Coef = e(b)
	mat Iter_Results = Coef[1,1],Coef[1,2],Coef[1,`=colsof(Coef)']
	mat rownames Iter_Results = `k'

	mat Results = Results\Iter_Results
}

// Save regression results for simulation
clear
set obs `iterations'
svmat Results, names(col)
save "$DATA/Temp/bootstrap_coefficients.dta", replace
*/


// 3.2 // Run simulation for each iteration
/*
use "$DATA/Temp/bootstrap_coefficients.dta", clear
mkmat b_1 b_2 constant, mat(Results)

use "$DATA/Temp/full_compressed.dta", clear
set more off

collapse (mean) rate gas p_spot, by(facility_id)

egen total_gas = sum(gas) // Total gas in sample
replace gas = gas*(32592000000/total_gas) // Scaling up gas production based on EIA 2016 estimate
gen ch4_new = rate*gas // Total emissions for calculating percentages later
egen total_ch4 = sum(ch4_new)
sum total_ch4
scalar total_emissions = r(mean)
drop ch4_new total_ch4 total_gas

egen lower_bound = min(rate)
 
order lower_bound, after(facility_id)
rename rate rate_0
rename gas gas_0
rename p_spot price_0
gen emissions_0 = rate_0*gas_0
gen ag_cost_0 = 0
gen ag_reduction_0 = 0
gen ag_value_0 = 0

gen rate_00 = rate_0 // Reference so others can be deleted
gen price_00 = price_0 // Reference so others can be deleted


compress *
save "$DATA/Temp/sim_base_compressed.dta", replace

local iterations $iterations
local sim_iterations 570 // Only running simulation up to just above social cost of methane for table results
local step_size .05

forval k = 1/`=`iterations'+1' {
	di `k' - 1
	use "$DATA/Temp/sim_base_compressed.dta", clear
	forval i = 1/`sim_iterations' {
		gen price_`i' = price_00 + .79*`step_size'*`i'
		gen rate_`i'= rate_`=(`i'-1)' + .79*`step_size'*(`=fp_pwr_1'*Results[`k',1]*price_`i'^`=fp_pwr_1-1' + Results[`k',2]*(`=fp_pwr_2'*(price_`i'^`=fp_pwr_2-1')*log(price_`i') + (price_`i'^`=fp_pwr_2')*(1/price_`i')  ) ) // First derivative
		
		qui replace rate_`i' = lower_bound if rate_`i' < lower_bound
		qui replace rate_`i' = rate_00 if rate_00 < lower_bound // Reference start point
		
		gen emissions_`i' = rate_`i'*gas_0
		gen reduction_`i' = emissions_`=(`i'-1)' - emissions_`i'
		gen value_`i' = reduction_`i'*price_00 // Value of conserved gas
		gen cost_`i' = ((price_`i' + price_`=(`i'-1)')/2)*(emissions_`=(`i'-1)' - emissions_`i')
		
		egen sum_cost_`i' = sum(cost_`i')
		egen sum_reduction_`i' = sum(reduction_`i')
		egen sum_value_`i' = sum(value_`i')
		gen ag_cost_`i' = ag_cost_`=`i'-1' + sum_cost_`i'
		gen ag_reduction_`i' = ag_reduction_`=`i'-1' + sum_reduction_`i'
		gen ag_value_`i' = ag_value_`=`i'-1' + sum_value_`i'
		gen mc_`i' = sum_cost_`i'/sum_reduction_`i'
		
		drop emissions_`=(`i'-1)' reduction_`i' value_`i' cost_`i' sum_* rate_`=(`i'-1)' price_`=(`i'-1)' //Dropping unnecessary variables to increase speed
		//graph addplot scatter rate_`i' price_`i' if rate_`i' < .02, msize(vsmall)
	}
	
	keep ag* mc*
	
	drop ag_cost_1-mc_60 ag_cost_70-mc_250 ag_cost_260-mc_560 // Only keeping the ones right around the values of interest to speed the simulation
	qui keep if _n == 1
	
	qui gen dummy = 1
	qui reshape long ag_cost_ ag_reduction_ ag_value_ mc_, i(dummy) j(tax)
	qui drop dummy
	
	rename tax tax_mcf
	rename mc_ mc_mcf
	rename ag_reduction_ ag_reduction_mcf
	rename ag_cost_ ag_cost
	rename ag_value_ ag_value
	
	qui replace tax_mcf = tax_mcf*`step_size'

	qui recast double tax_mcf
	qui replace tax_mcf = round(tax_mcf,.01)
	qui gen tax_tco2 = (tax*53.68)/34
	qui gen ag_reduction_tco2e = ag_reduction_mcf*34/53.68
	qui gen mc_tco2e = mc_mcf*53.68/34
	qui gen reduction_percent = (ag_reduction_mcf/total_emissions)*100
	qui gen net_cost = (ag_cost - ag_value)/28400000000
	
	qui replace ag_reduction_mcf = 0 if missing(ag_reduction_mcf)
	qui replace ag_reduction_tco2e = 0 if missing(ag_reduction_tco2e)

	format mc* tax* reduction %6.2fc
	format ag* %12.0fc

	order tax_mcf tax_tco2 mc_mcf mc_tco2e ag_cost ag_reduction_mcf ag_reduction_tco2e reduction_percent

	// Interpolating $5, $20, and $42 carbon taxes and social cost of methane
	qui set obs `=_N+1'
	qui replace tax_tco2 = 5 if missing(tax_tco2)
	qui set obs `=_N+1'
	qui replace tax_tco2 = 20 if missing(tax_tco2)
	qui set obs `=_N+1'
	qui replace tax_tco2 = 44.4 if missing(tax_tco2) // Social cost of methane
	sort tax_tco2
	
	qui gen percent_prev = (tax_tco2[_n+1] - tax_tco2)/(tax_tco2[_n+1] - tax_tco2[_n-1]) if missing(mc_mcf)
	qui gen percent_next = (tax_tco2 - tax_tco2[_n-1])/(tax_tco2[_n+1] - tax_tco2[_n-1]) if missing(mc_mcf)
	qui replace tax_mcf = tax_mcf[_n-1]*percent_prev + tax_mcf[_n+1]*percent_next if missing(tax_mcf)
	foreach var in mc_mcf mc_tco2e ag_cost ag_reduction_mcf ag_reduction_tco2e reduction_percent ag_value net_cost {
		qui replace `var' = `var'[_n-1]*percent_prev + `var'[_n+1]*percent_next if missing(`var')
	}
	
	qui keep if !missing(percent_next)
	qui gen iteration = `k' - 1
	qui drop percent_next percent_prev
	order iteration
	qui compress*
	
	qui save "$DATA/Temp/boot_stats_`=`k'-1'.dta", replace
}

*/


// 4 // Combine results
/*
use "$DATA/Temp/boot_stats_0.dta", clear

local iterations $iterations
drop tax_tco2 iteration
foreach var of varlist mc_mcf-net_cost {
	rename `var' `var'_0
}

save "$DATA/Temp/boot_stats_00.dta", replace

use "$DATA/Temp/boot_stats_0.dta", clear
forval k = 1/`iterations' {
	append using "$DATA/Temp/boot_stats_`k'.dta"
}


save "$DATA/Temp/boot_results_main.dta", replace

*/


// Construct outputs
/*
use "$DATA/Temp/boot_results_main.dta", clear

sort iteration tax_tco2 

gen avg_cost_tco2e = (ag_cost - ag_value)/ag_reduction_tco2e
format avg_cost_tco2e %6.2fc

recast double tax_tco2
replace tax_tco2 = 44.4 if tax_tco2 > 44.39 & tax_tco2 < 44.41 // Prevoiusly 44.400002 due to a variable formatting bug

gen double ag_reduction_tco2e_round = round(ag_reduction_tco2e,1000)
format ag_reduction_tco2e_round %15.0fc
drop ag_reduction_tco2e
rename ag_reduction_tco2e_round ag_reduction_tco2e

gen total_net_cost = (ag_cost - ag_value)/1000000
gen climate_benefit = (ag_reduction_tco2e*44.4)/1000000
gen welfare = climate_benefit - total_net_cost

gen emitted_gas = ag_reduction_tco2e*(1/(reduction_percent*.01))*(1-reduction_percent*.01)
gen taxes_collected = emitted_gas*tax_tco2
gen tax_per_mcf = taxes_collected/28400000000
gen cost_to_firms_per_mcf = tax_per_mcf + net_cost
format climate_benefit emitted_gas taxes_collected welfare %14.0fc
format tax_per cost_to %6.4fc

replace ag_cost = ag_cost/1000000
replace ag_value = ag_value/1000000

file open sim_results_wf using "$DATA\Results\sim_results_welfare.txt", write replace

foreach tax in 5 20 44.4 {
	sum tax_mcf if _n <= 3 & tax_tco2 == `tax' // Mean is iteration 0
	file write sim_results_wf %4.2fc (r(mean)) " & "
	
	sum tax_tco2 if _n <= 3 & tax_tco2 == `tax'
	file write sim_results_wf %4.2fc (r(mean)) " & "
	
	sum reduction_percent if _n <= 3 & tax_tco2 == `tax'
	file write sim_results_wf %4.1fc (r(mean)) "   & "

	sum ag_reduction_tco2e if _n <= 3 & tax_tco2 == `tax'
	file write sim_results_wf %10.0fc (r(mean)) "   & "
	
	sum ag_cost if _n <= 3 & tax_tco2 == `tax'
	file write sim_results_wf %3.0fc (r(mean)) "  & "

	sum ag_value if _n <= 3 & tax_tco2 == `tax'
	file write sim_results_wf %3.0fc (r(mean)) "  & "
	
	sum total_net_cost if _n <= 3 & tax_tco2 == `tax'
	file write sim_results_wf %3.0fc (r(mean)) "  & "
		
	sum net_cost if _n <= 3 & tax_tco2 == `tax'
	file write sim_results_wf %6.4fc (r(mean)) "   & "
	
	sum climate_benefit if _n <= 3 & tax_tco2 == `tax'
	file write sim_results_wf %5.0fc (r(mean)) "   & "
	
	sum welfare if _n <= 3 & tax_tco2 == `tax'
	file write sim_results_wf %5.0fc (r(mean)) "  \\" _n

	
	
	sum reduction_percent if _n > 3 & tax_tco2 == `tax' // Bootstrapped SEs
	file write sim_results_wf "     &      & (" %4.1fc (r(sd)) ") & ("	
	
	sum ag_reduction_tco2e if _n > 3 & tax_tco2 == `tax'
	file write sim_results_wf %10.0fc (r(sd)) ") & ("	
	
	sum ag_cost if _n > 3 & tax_tco2 == `tax' 
	file write sim_results_wf %2.0fc (r(sd)) ") & ("	

	sum ag_value if _n > 3 & tax_tco2 == `tax' 
	file write sim_results_wf %2.0fc (r(sd)) ") & ("
	
	sum total_net_cost if _n > 3 & tax_tco2 == `tax'
	file write sim_results_wf %2.0fc (r(sd)) ") & ("
		
	sum net_cost if _n > 3 & tax_tco2 == `tax'
	file write sim_results_wf %6.4fc (r(sd)) ") & ("
	
	sum climate_benefit if _n > 3 & tax_tco2 == `tax'
	file write sim_results_wf %3.0fc (r(sd)) ") & ("
	
	sum welfare if _n > 3 & tax_tco2 == `tax'
	file write sim_results_wf %3.0fc (r(sd)) ") \\" _n "[.75em]" _n
}
file close sim_results_wf

